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Strategies  for  automatic  design  of  power-optimized  3D-microbattery  geometries  are  here  investigated  by 
utilization  of  the  level-set  method  with  structure  topology  optimization.  The  methodology  is  extended 
from  solid  mechanics  to  electrochemical  systems,  where  battery  operation  is  simulated  using  the  Nernst 
—Planck  equation.  The  calculations  are  carried  out  for  the  3D-"trench”  geometry  with  LiCoCh  and  LiC6  as 
electrodes,  separated  with  a  LiPFe-PECho  polyethylene  oxide  polymer  electrolyte.  With  the  goal  to  ach¬ 
ieve  a  maximum  uniform  electrochemical  activity  over  the  electrode  surface  area,  an  optimized  electrode 
design  is  produced  by  coating  the  current  collectors  non-uniformly  with  active  material.  This  is  shown  to 
be  an  effect  of  the  3D  design  of  the  cell.  Evaluation  of  the  resulting  optimized  cell  by  simulations  of  the 
discharge  process  demonstrates  uniform  electrode  material  utilization  and  almost  uniform  current 
density  distribution  over  the  entire  electrode— electrolyte  interface.  Comparisons  between  optimized  and 
non-optimized  geometries  showed  that  the  geometry  optimization  increased  the  cell  performance  up  to 
2.25  times.  This  effect  is  mainly  achieved  by  minimizing  the  internal  energy  losses  caused  by  non¬ 
uniformities  in  the  ionic  transport  in  the  battery. 

©  2012  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Several  emerging  microtechnology  branches,  such  as  micro¬ 
electromechanical  system  (MEMS)  devices  and  small-scale  medical 
devices  like  e-pills,  require  portable  power  supplies  of  ~1  mm3 
dimensions  or  smaller.  Unfortunately,  the  development  of  such 
miniature  portable  power  sources  has  been  outpaced  by  the 
microtechnology  development,  resulting  in  technical  devices 
which  lack  appropriate  power  storage  [1],  To  fill  this  gap,  small 
scale  high  power  and  energy  density  batteries,  capable  of  providing 
both  high  current  and  capacity,  are  needed  [1],  Therefore,  a  new 
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type  of  thin-film  microbatteries  with  3D  electrodes  —  the  3D- 
microbatteries  (3D-MBs)  -  has  been  suggested  [1—3].  The  required 
combination  of  high  power  and  energy  density  is  achieved  by 
constructing  the  electrodes  in  3D-arrangements,  generating  large 
internal  surface  areas,  high  electrode  capacities  and  small  battery 
footprint  areas. 

The  manufacturing  of  these  kinds  of  electrodes  is  technically 
complicated,  thereby  often  resulting  in  electrodes  with  non- 
uniform  material  coating  [3].  At  the  same  time,  the  very  first 
mathematical  simulations  of  the  3D-MB  systems  have  demon¬ 
strated  non-uniform  current  densities  in  the  electrodes  and  elec¬ 
trolytes  [4-6],  indicating  possible  performance  losses  due  to 
inhomogeneous  electrode  reactions.  It  is  then  likely  that  these  two 
different  non-uniformities  co-interact  in  the  battery.  Previous 
simulations  [6]  have  demonstrated  how  uniformly  coated  elec¬ 
trodes  suffer  performance  losses  since  certain  parts  of  the  material 
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are  discharged  or  charged  faster  than  others.  This  could,  to  some 
degree,  be  adjusted,  for  example,  by  using  different  electrolytes  [7], 

An  optimal  3D-MB  performance  with  respect  to  both  power  and 
energy  supply  is  achieved  when  the  electrochemical  activity  over 
the  entire  electrode  surfaces  is  uniform  during  battery  operation. 
One  way  to  obtain  such  a  situation  would  be  to  utilize  electrodes 
with  a  non-uniform  material  layer  thickness  distribution  and  fit  it 
to  the  non-uniformity  in  the  current  distribution  resulting  from  the 
3D-design.  Considering  the  complexity  of  the  system,  advanced 
optimization  methods  are  needed  for  such  calculations. 

Structural  optimization  methods  with  compliance  minimization 
from  structural  mechanics  [8]  present  a  suitable  approach,  since 
these  methods  allow  a  systematical  and  simultaneous  design  of  the 
electrode  geometry  toward  its  optimal  current  distribution. 
However,  homogenization  methods  such  as  solid  isotropic  material 
with  penalization  (SIMP)  [9]  can  create  geometries  with  large  voids 
which  are  unsuitable  for  electrode  material  deposition.  A  region- 
based  approach,  such  as  the  level  set  method  [10—12]  is  therefore 
preferable. 

The  level  set  method  is  often  used  to  solve  compliance  mini¬ 
mization  problems  in  structural  mechanics  [11—13],  or  to  model 
multiphase  flow  in  fluid  dynamics  [14,15],  Recently,  a  substantial 
number  of  studies  applying  the  level  set  method  for  physical 
problems  have  been  published.  Zhuang  et  al.  [16]  used  the  level  set 
method  to  optimize  heat  transport  problems,  while  Luo  et  al.  [17] 
optimized  the  shape  and  topology  of  an  electromechanical 
actuator. 

In  the  current  study,  the  modeling  of  the  3D-MB  in  the  “trench” 
geometry  [1,3,6]  has  been  carried  out  using  the  Nernst— Planck 
equation  [18]  and  concentrated  solution  theory  [19—24],  The 
simulations  serve  two  major  goals  —  to  understand  the  ionic 
transport  in  the  3D  battery  system  and  to  optimize  the  3D-MB 
geometry.  Therefore,  an  optimization  algorithm  utilizing  the  level- 
set  method  is  proposed  for  a  3D-MB,  and  a  suitable  objective 
function  for  optimization  is  provided.  The  3D-MB  geometry  is 
thereafter  optimized,  and  simulations  of  discharge  processes  are 
carried  out  to  compare  the  optimized  and  non-optimized  cells. 


2.  Materials  and  methods 

2.1.  The  level  set  method 

Generally  in  the  level  set  method,  a  moving  material  boundary 
is  expressed  as  a  zero  level  set  of  a  higher  dimensional  function  in 
the  design  domain.  This  is  schematically  illustrated  in  Fig.  1.  The 
level  set  function  < p  divides  the  material  region  into  two  different 


materials,  so  that: 

(j>(x)  >  0,  x&Q\  (region  1)  (1.1) 

(j>(x)  =  0,  xedQ  (region  boundary)  (1.2) 

(/>(x)  <  0,  xeQ2  (region  2)  (1.3) 

Then,  different  regions  in  the  design  domain  are  defined  by  the 
Heaviside  function: 

m-{l:  (2) 

For  example,  the  spatial  conductivity  distribution  of  the 
electrode-electrolyte  system  is  expressed  as: 

o(<t>)  =  o,-H{4,)  +  o2-{1-H{4>))..  (3) 


function.  The  zero  level  set  represents  the  material  boundary  and  divides  the  whole 
region  into  two  subregions  with  different  material  properties.  When  the  level  set 
function  is  modified  during  the  calculation,  the  material  boundary  is  moved. 


where  a\  represents  the  (electronic)  conductivity  of  the  electrodes 
and  <72  is  the  (ionic)  conductivity  of  the  electrolyte.  In  conjunction, 
these  describe  the  total  conductivity  distribution  in  the  cell. 

To  avoid  numerical  problems  at  the  region  interfaces, 
a  smoothed  Heaviside  function  is  used  during  the  calculations 
[11,13],  The  boundary  of  the  regions  is  identified  by  using  the  Dirac 
delta  function,  defined  as  <5 (</>)  =  and  expressed  depending 
from  the  level  set  function: 

(4) 

The  normal  vector  at  the  region  interfaces  is  defined  as: 
n  =  V0/|V0|,  (5) 

Finally,  the  level  set  function  is  a  solution  of  Hamilton  Jacoby 
type  equation: 

Vn|V0|  =  0,  (6) 

where  Vn  is  the  normal  velocity  of  the  moving  boundaries  in  the 
system.  By  solving  this  equation,  it  is  possible  to  track  the  move¬ 
ment  of  the  regions  in  the  design  domain. 

2.2.  The  objective  function 

The  strategy  to  optimize  the  3D-MB  geometry  is  to  reach  as 
uniform  electrochemical  activity  as  possible  on  the  electrode 
surface.  While  this  condition  is  fulfilled,  the  lithiation/delithiation 
of  the  electrodes  is  uniform  and  the  battery  performance  is 
maximal. 
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To  carry  out  an  optimization  procedure,  it  is  necessary  to  define 
a  function  to  be  minimized  (or  maximized)  —  an  objective  function. 
Here  we  consider  an  objective  function  based  on  the  current 
density  in  the  electrolyte,  which  accounts  for  the  main  resistance  in 
the  battery  cell.  The  current  density  is  directly  correlated  to  the 
concentration  dependent  ionic  conductivity,  and  is  a  sum  of  the 
diffusive  and  migrative  fluxes. 

Similarly  to  the  optimization  of  heat  transport  [16],  we  propose 
an  objective  function  as  a  functional  of  the  current  density, 
expressed  as: 

F=  f  +Py  +J?d  Q  (7) 

The  objective  function  defined  in  Eq.  (7)  will  be  minimized,  if 
the  average  current  density  in  the  cell  is  as  uniform  as  possible. 


value.  The  solution  of  the  optimization  problem  is  terminated  if  the 
following  condition  is  fulfilled  [11]: 

|[Vn|^)|V0|df3<rtol  (11) 

where  ytoi  is  specified  error  limit. 

To  calculate  the  Lagrange  multiplier,  the  volume  of  the  elec¬ 
trodes  must  be  constant  [11]: 


(12) 

Now,  the  Lagrange  multiplier  can  be  expressed  as  [11]: 


2.3.  General  optimization  problem 

By  minimizing  the  objective  function,  the  most  uniform  current 
density  distribution  in  the  3D-MB  geometry  is  achieved.  To  ensure 
that  the  battery  electrodes  will  not  fill  the  whole  cell,  a  constraint  is 
introduced:  the  battery  electrodes  sustain  a  specific  volume  V*.  The 
constraint  is  introduced  by  using  the  Heaviside  function  applied  to 
the  level-set  variable.  To  minimize  the  objective  function,  specified 
in  Eq.  (7),  the  general  topology  optimization  problem  for  the  3D-MB 
can  then  be  written  as  [8]: 

minimize  F(</>)  =  J  +  jj  +  jjdQ  (8.1) 

Q 

subject  to  J  H(0)dfi  =  V*  (8.2) 


where  V*  is  the  volume  of  the  electrodes  and  F  is  the  objective 
function,  calculated  by  a  system  of  partial  differential  equations 
presented  in  paragraph  2.5. 

The  implementation  of  the  optimization  routine  follows  Wang 
et  al.  [11],  Accordingly,  the  constraints  in  the  optimization  are 
implemented  using  the  method  of  Lagrange  multipliers: 

Wb  )+<HWiOdfi  (9) 

The  speed  function  is  achieved  by  minimizing  Eq.  (9),  as 
demonstrated  by  Wang  et  al.  [11].  Since  the  mathematical  deriva¬ 
tion  of  the  final  speed  function  can  be  found  in  the  literature  [11] 
and  it  is  not  an  aim  of  current  work,  we  present  here  only  the 
final  result.  The  normal  velocity  of  the  zero  level  set  of  the  level  set 
surface  can  be  calculated  as  Vn  =  (F+  A)<5(0)  [11,13].  It  is  known  from 
structural  mechanics  that  it  is  sufficient  to  know  the  normal 
velocity  of  the  domain  boundary  to  describe  the  change  of  its 
shape.  As  a  result,  by  using  Eq.  (6),  the  evolution  of  the  electrode 
boundary  toward  its  optimal  shape  can  be  calculated  [11,13]: 

d^  +  (F  +  X)8(m\  =  0,  (10) 

where  t  is  pseudo  time,  representing  the  optimization  steps. 
Pseudo  time  is  used  instead  of  time  to  emphasize  that  while  the 
shape  of  the  electrodes  evolves  toward  optimal,  the  mass  and 
current  distribution  in  the  cell  are  at  steady  state  during  every 
optimization  step.  When  solving  this  optimization  problem,  the 
solution  will  converge  toward  a  local  minimum  based  on  the  initial 


J  F<52(</>)|</>|dfl 

A  =  -  (13) 

j  t>2(m\M 


2.4.  Solution  of  the  level  set  equation 

The  level  set  equation  is  solved  in  Comsol  Multiphysics  4.2  [25] 
by  using  the  “Level  Set  toolbox”  in  the  computational  fluid 
dynamics  module.  The  toolbox  includes  methods  for  reinitializa¬ 
tion  and  numerical  stability:  the  reinitialization  is  needed  to  keep 
the  level  set  function  a  signed  distance  function  with  the  property 
that  the  norm  of  the  gradient  of  the  level  set  function  is  as  close  to  1 
as  possible  (within  numerical  errors)  [26],  In  current  work,  the 
reinitialization  is  introduced  according  to  Olsson  et  al.  [14,15]  since 
the  Level  Set  toolbox  in  Comsol  Multiphysics  is  used.  During  the 
solution,  the  region  boundaries  are  represented  by  the  <p  =  0.5  level 
set  value.  The  equation  provided  by  Olsson  et  al.  [14,15]  describes 
the  evolution  of  the  level  set  field  depending  on  the  velocity  vector 


U-V0  =  TV-(£lsV0- 0(1-0^)  (14) 

Here,  u  —  u{x,y )  is  the  velocity  of  the  level-set  surface.  To  solve 
the  optimization  problem,  the  normal  velocity  is  needed.  By  using 
the  normal  vector  n  and  normal  velocity  l/n  of  the  surface,  the 
velocity  u  can  be  expressed  as  it  =  Vn-3t  =  Vn  •  (v$/|  V0|). 

As  a  result  the  "u  V0term  can  be  replaced  with  the  relation: 


leading  to  the  following  equation  to  be  solved: 


(15) 


^+Vn|V0|  =  yV-(qsV0-<Kl  0)  vjj)-  H6) 

In  equation  (16),  the  left  hand  side  describes  the  movement  of 
the  level  set  surface  while  the  right  hand  side  guarantees  numerical 
stability. 


2.5.  Battery  materials  and  model 

In  the  current  study,  a  LiPFe-PEC^o  polyethylene  oxide  electro¬ 
lyte  equivalent  to  1.5  M  LiPF6  is  simulated.  The  battery  electrode 
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materials  are  UC0O2  (positive  electrode)  and  graphite  (negative 
electrode). 

In  the  calculations,  a  3D-MB  trench  geometry  is  simulated 
(Fig.  2a)  consisting  of  planar  Cu  current  collectors  bases  with 
attached  perpendicular  A1  current  collector  plates,  both  coated  with 
electrode  materials.  The  trench  geometry  it  is  one  of  the  most 
promising  3D-MB  geometries  [3]  and  can  be  considered  a  first  step 
toward  studying  geometrically  more  complex  3D-geometries.  In 
the  current  work,  the  battery  simulation  follows  well  established 
and  experimentally  validated  electrochemical  descriptions  [18— 
24],  Material  properties  such  as  diffusion  coefficients  and  conduc¬ 
tivities  of  the  electrodes  and  electrolytes  are  expressed  by  the  level- 
set  function,  making  it  possible  to  describe  changes  in  the  electrode 
shape  during  the  optimization.  During  the  optimization  process, 
the  battery  is  simulated  at  steady-state,  so  that  the  material  balance 
in  the  electrodes  does  not  need  to  be  taken  into  account.  The 
Nernst-Planck  equation  is  used  to  model  the  mass  and  current 
transport  in  the  electrolyte,  since  previous  studies  [18]  have 
demonstrated  excellent  agreement  between  experimental  and 
simulation  results  at  constant  current.  Furthermore,  the  Nernst- 
Planck  equation  is  more  convenient  to  use  than  concentrated 
solution  theory  based  models  [19-24]  during  the  optimization 
routine  described  above,  since  it  is  computationally  less  expensive 
and  more  robust  [27], 

When  the  material  boundary,  expressed  by  the  zero  level  set 
of  the  level  set  function,  moves  through  the  mesh  elements,  it  is 
then  possible  to  use  only  homogeneous  Neumann  boundary 
conditions  or  continuity  boundary  conditions  for  both  electrical 
current  and  ionic  flux  at  the  electrode— electrolyte  interface,  since 
the  material  boundary  is  generally  not  conforming  with  the  mesh 
element  edges.  In  order  to  use  previously  validated  battery 
models  [18-24],  porous  electrodes  must  be  implemented,  since 
homogeneous  Neumann  boundary  conditions  then  can  be  applied 
at  the  electrode— electrolyte  interfaces.  The  use  of  porous  elec¬ 
trodes  is  also  advantageous  since  random  numerical  errors  may 
occur  when  the  material  boundaries  move  through  the  elements, 
which  can  stop  the  optimization  iterations  prematurely  [13]. 
These  random  errors  are  smoothed  out  when  utilizing  porous 
electrodes. 


To  perform  the  calculations,  the  diffusion  coefficients  and 
conductivities  must  be  expressed  through  the  level  set  function. 
Consequently,  the  conductivity  of  the  electrolyte  is 

<h(0)  =  1  (17) 

where  the  first  term  in  the  right  side  represents  the  conductivity  of 
electrolyte  inside  the  porous  electrode  material  and  the  second 
term  the  conductivity  of  pure  electrolyte.  The  conductivity  in  the 
electrodes  is 

=  (l-A'HW  +  ^-a-^))-  J  =  n>  P  (18) 

where  j  =  n  corresponds  to  positive  electrode  and  j  —  p  corresponds 
to  negative  electrode,  ao  is  an  artificial  parameter  of  very  small 
value,  introduced  to  avoid  numerical  singularities  [9],  Finally,  the 
diffusion  coefficients  in  the  electrolyte  are  formulated  by: 


Dj($)  =  El5DiH(4>)  +D,-(1  -H(0)),  |  =  Li,  PF6,  (19) 

The  Butler-Volmer  equation  describes  the  charge  transport 
between  the  electrode  and  electrolyte: 


where  p  =  cps  -  <p\  -  Voc.  The  open  circuit  potential  is  fitted  to 
experimental  data  according  to  Ref.  [28],  Since  the  electrode  reac¬ 
tions  in  Li  batteries  are  fast  and  usually  limited  by  mass  transport, 
a  high  constant  exchange  current  density  (i'o  =  100  A  m-2)  is  used 
during  the  optimization  to  avoid  numerical  errors  which  may  arise 
at  the  electrode— electrolyte  boundary,  expressed  by  the  level-set 
function.  During  the  time  dependent  simulations  of  the  discharge 
processes,  the  exchange  current  density  is  calculated  as 
i0  =  F(/ca)ac(/Oac(4m^  where  j  =  n,p, 

k  =  ka  =  kc  =  2  •  10_1*  m  s_1  and  aa  =  ac  =  0.5  in  order  to  achieve  an 
initial  exchange  current  density  close  to  values  reported  in  litera¬ 
ture  [23,24], 


base  Al  ^ 


Fig.  2.  Geometry  c 
the  test  problem. 


:  3D-trench  i 


I  electrode  plate  arrangement  in  3D-trench  geometry;  (a)  represents  the  geometry  of  the  optimization  problem  and  (b)  the  geometry  of 
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When  optimizing  the  electrode  shape,  the  flux  of  specie  i  obeys 
the  Nernst— Planck  equation  [29]: 


Ji  =  -AWVq-^D.WCjV^  i  =  Li,  PF6  (21) 

where  q  is  the  concentration,  D,-  the  diffusion  coefficients  of  species 
i,  Zj  the  charge  number,  F  Faraday’s  constant,  R  the  universal  gas 
constant  and  <p  the  potential  in  electrolyte.  Electro-neutrality  in  the 
electrolyte  is  assumed,  so  that  c  =  cu  =  cPFf.. 

The  concentration  profile  at  steady  state  can  then  described 
according  to  Ref.  [18],  however,  since  the  porous  electrodes  are 
used,  the  Li  insertion  is  described  by  using  the  source  term  [24,30]: 


2DLiWDPF6W  2 
(Du(0)+Dpf6(0)) 


l-t° 

¥ 


ad 


(22) 


The  initial  value  of  the  concentration  (Co)  is  specified  in  Table  1. 
Since  Li+  ions  can  be  inserted  or  removed  from  the  electrolyte  only 
by  the  electrode  reactions,  the  boundary  condition  for  Eq.  (10)  is 
[18]: 


2DLiWDpF6W  ^  = 

(Du(0)+Dpf6(0)) 


(23) 


The  steady-state  solution  of  Eq.  (10)  requires  an  extra  constraint 
to  guarantee  a  unique  solution.  This  is  achieved  by  specifying  the 
conservation  of  mass  as  [7]: 


Table  1 

Parameters  used  in  the  simulations. 


Symbol  Quantity 

Du  Diffusion  constant  for  Li+  ions  in  polymer 
electrolyte  [33] 

Dpp,  Diffusion  constant  for  PFg  ions  in  polymer 
electrolyte  [33] 

D  Li  diffusion  constant  in  the  active  material  [34] 

a ai  Electronic  conductivity  of  positive  current 

collector  (Al)  [25] 

o-p  Electronic  conductivity  of  positive  electrode 

active  material 
(UCo02)  [34] 

ffn  Electronic  conductivity  of  negative  electrode 

(UC6)  [34] 

<rCu  Electronic  conductivity  of  negative  current 
collector  (Cu)  [25] 

c0  Salt  concentration  in  electrolyte  [7] 

jo  Discharging  current  during  optimization 

r  Porosity  of  the  electrodes 

i0  Exchange  current  density  during  optimization 

k  Rate  coefficient  during  discharging  simulations 

positive  electrode 

c"max  Maximum  Li  concentration  in  the  active 
material  of  the 
negative  electrode 
H  Height  of  the  plate 

h  Distance  between  plate  tip  an  opposing 

electrode  base 

d  Distance  between  current  collector  plates 

dc  Thickness  of  the  level-set  function  separating 

%  Parameter  controlling  electrode-electrolyte 

interface 

thickness  during  optimization  (during  validation) 
y  Reinitialization  parameter  during  optimization 

(during  validation). 


2.5  10-13m2s-' 
3  10~13  m2  s-1 

2.5 -10-13  m2  s"1 
3.75  107  S  m_1 

MO^Sm-1 


ISm-' 

5.95 -107  S  ITT1 

1.5  M 
10  A  nr2 
0.5 

100  A  nr2 
2-10~n  m  s_1 
51,656  mol  m  3 


28,225  mol  m-3 


40- 10  6 m 
15  10~6  m 

15106  m 
110  6m 

0.5  10~6  m 
(1  10~6  m) 

100  m  s-1 
(40  m  s1) 


Here,  12  is  the  electrolyte  region,  V  the  volume  of  the  electrolyte 
and  Co  the  initial  concentration. 

During  the  optimization,  the  state  of  charge  is  kept  constant. 
Moreover,  the  Li  concentration  in  the  electrodes  is  uniform,  and  it  is 
assumed  that  the  positive  electrode  is  almost  completely  lithiated 
(cs  =  cPmax  -  1000  mol  m-3)  and  negative  almost  completely 
delithiated  (cs  =  1000  mol  m-3),  corresponding  to  a  discharged  cell. 
This  discharged  cell  is  used  since  the  result  will  guarantee  uniform 
material  utilization. 

The  potential  in  the  electrolyte  is  calculated  according  to  Refs. 
[7,27,29],  since  the  porous  electrodes  are  used,  reaction  current  is 
specified  by  the  source  term  [24,30]: 

(Dpf6(<«  -  du(0))fv2c  -  V-  [  J  (Dpp6(0)  +  Du(fl) V*]  =  ad 

(25) 


n>-[(DpF6(0)-Du(0))FVc-^(DpF6W+DLi(0))V¥,]  =  0 

(26) 

where  as  is  the  specific  surface  area  of  the  electrodes  and  is 
calculated  as  as  =  3c/rp,  where  rp  is  radius  of  an  active  material 
particle.  The  current  densities  in  the  electrodes  and  in  the  current 
collectors  are  calculated  by: 


Vo-,(0)V</)  =  ad,  i  =  Al,p,  n,  Cu 

(27) 

rt>'(ffAf(0)V<P)  —  Jo 

(28) 

<Plnegative_current_collector  =  0, 

(29) 

where  <p  is  the  electrode  potential,  jo  is  the  discharging  current  in  the 
current  collector  of  the  positive  electrode  and  <7;  is  the  conductivity. 
The  values  of  coefficient  i  represent  the  positive  current  collector 
(i  =  Al),  the  positive  active  material  (i  =  p),  the  negative  active 
material  (i  =  n)  and  the  negative  current  collector  (i  =  Cu). 

2.6.  Cell  performance  evaluation 

To  evaluate  the  performance  of  the  cell,  the  dissipated  energy 
during  the  discharge  process  is  calculated.  First,  the  overpotential 
of  the  cell  is  defined  as: 

AV7  =  Vdisc  -  Voc  (30) 

where  l/disc  is  the  discharging  voltage  and  Voc  is  the  open  circuit 
voltage  of  the  cell,  calculated  as  Voc  =  VuctiOj^p)  -  luCsO^n), 
where  xp  is  x  in  LixCo02  and  xn  is  x  in  LixC6  [28],  A  zero  overpotential 
represents  the  case  where  no  net  current  is  moving  through  the 
cell,  while  any  non-zero  overpotential  value  represents  the  energy 
required  to  operate  the  cell  at  a  specific  current. 

According  to  Bard  [29],  the  cell  resistance  can  be  calculated  as: 

r  =  AV/jo,  (31) 

where  J  is  the  discharge  current  density  and  r  is  the  combined 
charge  and  mass  transport  resistances.  The  use  of  current  density, 
instead  of  current,  results  in  a  normalized  resistance  r  (C2  m2). 
However,  a  numerical  value  of  the  resistance  can  be  calculated  by 
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multiplying  the  normalized  resistance  with  the  surface  area  of  the 
current  collector  (i.e.,  the  cell  fingerprint  area).  We  here  define  the 
cell  performance  (resulting  from  the  cell  geometry  optimization)  as 
the  ratio  of  dissipated  energy  during  discharge  in  the  reference  cell 
and  in  the  optimized  cell.  The  dissipated  power  is  calculated  by 
Joule— Lenz  law,  and  therefore  the  performance  increase  is: 

fnonopt  _  Jo'fnonopt  _  fnonopt  (32) 

fopt  Jq  *  t°pt  Popt 

where  rnonopt  is  the  resistance  in  the  reference  cell,  ropt  is  the  resis¬ 
tance  in  the  optimized  cell  and  jo  is  the  discharge  current  density.  The 
performance  increase  characterizes  the  increase  in  the  battery  power, 
since  the  volume  of  the  electrodes  (battery  capacity)  is  not  changing 
during  the  optimization  while  the  cell  resistance  is  decreased. 

2.7.  Simulation  setup 

2.7.1.  Validation  of  the  optimization  algorithm 

To  validate  the  optimization  routine,  a  suitable  test  problem  with 
known  optimal  solution  was  solved.  The  most  uniform  current 
density  in  an  electrochemical  cell  would  be  achieved  in  a  conventional 
2D-cell  (however,  with  limited  energy  storage  capability)  with  planar 
electrodes  having  equal  surface  area.  In  the  test  problem,  the  initial 
3D-MB  architecture  uses  flat  current  collectors,  as  illustrated  in  Fig.  2b. 
The  simulation  was  carried  out  in  the  area  between  the  dashed  lines, 
while  the  symmetry  of  the  structure  was  used  to  extend  the  result  to 
the  entire  battery  cell.  Since  the  electrode  and  electrolyte  spatial 
distribution  are  described  by  the  level-set  function  (Eqs.  (1)— (3)  and 
Fig.  1 )  and  the  level-set  function  value  1  corresponds  to  the  electrode, 
the  level  set  function  was  required  to  be  1  in  the  current  collectors  to 
ensure  that  at  least  some  amount  of  the  electrode  material  is  always 
attached  to  them  (Eq.  (18)).  The  optimization  algorithm  would  then 
converge  when  the  electrode  geometry  becomes  2-dimensional,  since 
a  homogeneous  current  density  is  achieved. 

In  the  calculations,  the  geometry  was  discretized  using  ~  15,000 
triangular  linear  elements. 

2.7.2.  Optimization  of  the  3D-trench  geometry  and  evaluation  of 
the  results 

The  optimization  of  the  electrode  shape  was  carried  out  using  the 
geometry  represented  in  Fig.  2a.  At  the  start  of  the  optimization 
procedure,  the  current  collectors  were  coated  with  a  uniform  layer  of 


electrode  materials  while  the  level  set  function  value  in  the  electrodes 
and  the  electrolyte  was  1  and  0,  respectively,  so  that  the  electrode- 
electrolyte  boundary  was  represented  by  a  0.5  value.  During  the 
optimization,  the  shape  of  both  electrodes  was  changed.  However, 
since  one  “optimal”  configurations  is  achieved  when  the  electrodes 
are  connected  (representing  a  short-circuited  cell),  a  thin  separating 
region  (dc  in  Fig.  2)  where  the  level  set  function  was  required  to  be 
0  was  introduced.  The  electrode  boundary  can  thereby  not  move  over 
this  region  during  the  optimization  process.  All  material  properties 
and  geometrical  parameters  are  specified  in  Table  1.  During  the 
optimization,  a  current  density  of  10  A  m-2  was  used. 

In  the  calculations,  a  triangular  mesh  with  ~  48,000  linear 
elements  was  used.  A  high  number  of  elements  were  required  to 
accurately  model  the  battery  regions  where  thin  electrode  material 
layers  were  developing. 

After  the  optimization  procedure  was  finished,  the  effects  of  the 
optimization  of  the  battery  geometry  was  evaluated  by  full  cell 
charge/discharge  simulations  and  compared  to  the  pre-optimized 
3D-trench  model.  Time  dependent  charge/discharge  simulations 
were  carried  out  with  the  optimized  geometry  and  compared  to 
simulations  of  a  reference  scenario.  The  reference  scenario  was  the 
starting  point  of  the  optimization  -  i.e.,  the  3D-trench  architecture 
with  uniform  electrode  thickness.  The  effect  of  the  optimization 
was  evaluated  by  discharging  fully  charged  cell;  the  current 
densities  used  for  evaluation  in  both  optimized  and  reference 
geometries  were  either  2  A  m-2,  4  A  itT2,  6  A  m-2,  8  A  m-2  or 
10  A  m-2.  For  the  evaluation  of  the  optimization  results,  new 
geometrical  model  of  the  optimized  battery  geometry  was  con¬ 
structed,  so  that  the  electrodes  were  modeled  geometrically, 
independently  from  the  level-set  function. 

In  the  calculations,  a  triangular  mesh  of  ~2700  linear  elements 
was  used  to  simulate  the  optimized  geometry.  For  the  reference 
geometry,  ~1700  triangular  quadratic  elements  were  used. 

For  the  time  dependent  charging— discharging  simulations,  the 
Battery  toolbox  in  Comsol  Multiphysics  4.2  [25],  utilizing  concen¬ 
trated  solution  theory  [19—24]  was  used. 


3.  Results  and  discussion 

3.1.  The  test  problem 

The  progression  of  the  solution  of  the  test  problem  described  in 
paragraph  2.7.1  is  presented  in  Fig.  3.  Black  color  represents  the 


optimization 


End  of  the 
optimization 


Fig.  3.  Electrode  shape  development  for  the  optimization  test  problem  during  the  calculation  (a-d).  Black  color  represents  electrode  material  and  white 
During  the  optimization,  the  3D-electrode  geometry  (a)  is  transformed  to  a  2D-electrode  (d). 


1  color  the  electrolyte. 
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electrodes  and  white  represents  the  electrolyte;  both  electrodes 
have  a  constant  volume  during  the  optimization.  It  is  clearly  seen 
how  the  optimization  routine  starts  to  remove  material  from  the 
electrode  plates  and  increase  the  thickness  of  the  electrode  base.  As 
expected,  the  plates  are  removed  during  the  course  of  the  optimi¬ 
zation  and  the  material  is  carried  over  to  the  base  of  the  electrode, 
forming  planar  electrodes  with  only  negligible  distortions.  The 
solution  of  the  test  problem  clearly  demonstrates  the  expected 
behavior  of  the  optimization  algorithm  for  this  multiphysics 
problem. 


3.2.  Optimization  of  the  3D-trench  geometry 

The  results  of  the  optimization  of  the  electrode  material 
distribution  in  the  3D-trench  geometry  are  presented  in  Fig.  4. 
Again,  black  represents  the  electrode  and  white  the  electrolyte 
area.  As  opposed  to  the  test  problem,  fixed  3D-structured  current 
collector  substrates  are  coated  with  electrode  materials  (see  Fig. 
2a),  thereby  reassuring  a  3D-MB  as  the  optimization  result.  The 
optimization  routine  progresses  toward  a  homogeneous  current 
density  (by  minimizing  the  objective  function).  At  the  beginning 
of  the  optimization  process,  the  current  density  is  very  high  at 
the  tips  of  the  electrodes  and  low  at  the  electrode  base.  To  obtain 
the  most  uniform  current  density,  the  optimization  routine  starts 
to  rearrange  the  active  material  so  that  the  electrode  areas  with 
high  current  density  gain  more  material,  resulting  in  increased 
local  surface  area  and  decreased  local  current  density  (due  to  that 
the  material  amount  increase  in  the  electrode  tips).  Alternatively, 
by  surrounding  certain  parts  of  the  electrodes  with  a  larger  mass 
of  electrolyte  higher  local  resistivity  values  are  achieved  (the 
amount  of  material  decreases  at  the  electrode  base).  This  process 
continues  gradually  during  the  whole  optimization  cycle  and 
finally  creates  electrodes  where  almost  all  active  material  has 
been  transferred  from  the  bottom  of  the  electrodes  to  the  tips. 


The  optimal  distribution  of  material  is  different  for  the  positive 
and  negative  electrode.  During  the  optimization,  most  material  on 
the  positive  electrode  is  assembled  at  the  tip  of  the  plate  and  only 
a  small  amount  is  left  at  the  base  of  the  electrode.  The  optimal 
shape  of  the  negative  electrode  is  more  complex.  At  the  beginning 
of  the  optimization,  the  evolution  of  both  electrode  shapes  is 
similar,  but  after  ~  50%  of  the  calculation  is  finished,  a  sharp  corner 
starts  to  develop  at  the  tip  of  the  electrode.  Moreover,  the  optimi¬ 
zation  routine  starts  to  move  material  from  the  center  of  the  tip  of 
the  negative  electrode  (Fig.  4c).  As  a  result,  the  local  surface  area  of 
the  electrode  tip  is  increasing  and  local  current  density  decreases. 

The  explanation  for  this  behavior  is  that  since  one  of  the  optimal 
electrode  configurations  is  achieved  when  the  positive  and  nega¬ 
tive  electrodes  are  connected  (which  represents  a  short-circuited 
cell),  a  thin  region  dc  (Fig.  2)  has  been  created  between  the  elec¬ 
trodes,  dividing  the  electrolyte  into  two  volumetrically  equal  parts. 
While  the  optimization  of  the  positive  electrode  is  unaffected  by 
this  constraint,  it  influences  the  shape  of  the  negative  electrode  at 
the  end  of  the  calculation  by  creating  a  negative  electrode  with 
sharp  comers  (Fig.  4d).  In  future  studies,  the  model  could  be 
improved  by  addition  of  moving  mesh  interfaces,  making  it 
possible  to  move  the  separating  region  away  from  the  surface  of  the 
optimized  electrode.  This  would  also  allow  optimization  of  elec¬ 
trodes  with  very  short  separating  distances. 

During  the  optimization,  both  electrodes  lose  ~10%  of  their 
volume.  This  loss  of  material  is  caused  by  the  Lagrange  multiplier, 
presented  in  Eqs.  (12)  and  (13)  since  the  total  volume  of  the  material 
is  not  mathematically  represented  when  using  the  derivative  (Eq. 
(12)).  The  Lagrange  multiplier  instead  enforces  zero  material 
change.  However,  small  numerical  errors  will  always  occur  during 
the  iterative  computer  simulations  which  gradually  can  allow  the 
system  to  deviate  from  its  initial  volume.  Such  a  small  violation  of 
the  volume  constraint  is  often  reported  in  literature  [12,13]  discus¬ 
sing  optimization  with  the  level  set  method.  However,  the  impact  on 
the  results  should  be  limited  since  it  is  easily  identified  as 
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Fig.  4.  Evolution  of  the  electrode  shape  during  the  optimization  procedure.  The  uniform  material  coat 
represents  the  porous  electrode  materials;  white  electrolyte. 
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a  systematic  error  in  all  comparisons  and  the  results,  comparing  the 
discharge  curves  and  performance  are  plotted  in  normalized  scale. 

Although  there  exists  a  wide  range  of  proposed  3D-MB  geom¬ 
etries,  the  “trench”  geometry  has  been  the  focus  of  the  current 
studies.  Other  3D-geometries  such  as  pillar  architectures  in  inter- 
digitated  or  concentric  arrangements  [1]  would  also  be  important 
to  study  as  good  candidates  for  3D-MBs.  However,  the  trench 
geometry  is  one  of  possible  3D-MB  geometries  closest  to  realiza¬ 
tion,  and  it  can  also  be  considered  as  a  model  for  the  interdigitated 
and  concentric  geometries.  For  example,  if  the  simulated  trench 
geometry  is  rotated,  a  concentric  geometry  is  obtained,  or  if  the 
electrodes  in  the  current  geometry  are  rotated,  an  interdigitated 
pillar  geometry  is  achieved.  Due  to  symmetry  considerations,  the 
current  density  distribution  in  the  interdigitated  design  would 
probably  be  somewhat  less  uniform  than  in  the  trench  geometry, 
since  the  current  density  distribution  will  change  also  around  the 
pillar  cross-sections  [4,6],  An  optimization  should  therefore  lead  to 
a  similar  electrode  material  distribution  as  presented  here,  but  with 
a  slight  material  layer  thickness  variation  around  the  pillar  cross- 
section. 


3.3.  Performance  of  the  optimized  cell 

The  current  density  distributions  in  the  optimized  and  reference 
(starting)  3D-MB  geometries  are  presented  in  Fig.  5  during  the 
discharge  of  the  cell  for  a  current  of  10  A  m-2.  The  current  density 
distribution  in  the  reference  geometry  demonstrates  a  significant 
development  during  the  course  of  the  discharge.  At  the  beginning 
of  the  process,  at  t  =  0  s,  the  current  density  is  maximal  at  the  tip  of 
the  positive  electrode  and  most  of  the  current  is  moving  between 
electrode  plates.  As  the  simulation  progresses  (at  t  =  350  s;  Fig.  5e), 
the  current  density  maximum  moves  to  the  negative  electrode, 
similarly  as  for  liquid  electrolytes  in  the  same  geometry  [6].  At  the 
same  time,  the  ionic  transport  pathways  (represented  by  current 
density  streamlines)  are  becoming  longer,  since  the  electrode 
material  at  the  tips  is  being  depleted  due  to  higher  activity  in  these 
regions.  Mass  and  charge  transport  continues  through  the  elec¬ 
trode  bases,  resulting  in  longer  transport  distances,  and  thereby 
higher  cell  resistance.  At  the  end  of  the  discharge  (at  t  =  700  s; 
Fig.  4f)  the  maximum  current  density  can  still  be  found  at  the 
negative  electrode,  but  now  the  current  moves  almost  completely 
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between  electrode  bases.  The  ionic  transport  pathways  and  cell 
resistance  are  significantly  increased  and  the  ions  must  travel  up  to 
60  pm  to  reach  the  opposite  electrode,  instead  of  10-20  pm  at  the 
beginning  of  the  discharge.  Thus,  the  energy  lost  due  to  the 
transport  of  the  ions  is  increased.  To  balance  this  effect  in  the  non- 
optimized  cell,  the  ionic  conductivity  of  the  electrolyte  has  to  be 
much  higher. 

The  optimized  geometry,  on  the  other  hand,  displays  much  less 
time  dependent  changes  in  the  current  density  distribution.  It  can 
be  seen  in  Fig.  5a— c  how  the  current  density  distribution  remains 
almost  unchanged  through  the  simulated  discharge  process  —  the 
changes  are  very  small  and  appear  at  the  end  of  the  discharge.  The 
most  significant  changes  appear  at  the  tip  of  the  negative  electrode, 
and  are  visible  in  Fig.  5b  and  c  where  a  significant  increase  in  the 
current  density  (as  compared  to  Fig.  5a)  is  seen.  This  increase  of 
current  density  is  caused  by  the  locally  limited  electrolyte  volume 
electrochemically  associated  with  a  large  electrode  surface  area.  On 
the  other  hand,  this  “hole”  at  the  tip  of  the  electrode,  constructed 
by  the  optimization  procedure,  makes  it  possible  to  discharge  this 
rather  large  part  of  the  electrode  uniformly. 

The  current  density  in  the  optimized  geometry  is  dominated  by 
the  electrode  plates,  which  is  expected,  considering  that  much 


material  from  the  bottom  of  the  electrodes  has  been  removed 
during  the  optimization  procedure.  However,  this  re-arrangement 
of  active  material  obviously  makes  it  possible  to  achieve  a  much 
more  uniform  current  density  over  the  whole  active  material 
region. 

The  effects  of  the  optimization  can  also  be  seen  in  Fig.  6,  where 
the  concentration  distribution  in  the  electrodes  is  presented.  It  can 
be  seen  that  the  Li  concentration  in  the  active  material  in  the 
optimized  electrodes  changes  uniformly  during  discharge  -  the 
concentration  distribution  in  almost  all  parts  of  the  electrode 
changes  with  equal  speed.  In  contrast,  the  reference  geometry 
shows  significant  changes  in  Li  concentration  in  the  electrodes;  the 
charging/discharging  starts  from  the  plate  tips  and  progresses 
through  the  whole  electrode  plate  toward  the  base,  which  is 
charged/discharged  last.  This  kind  of  charging/discharging 
dynamics  increases  the  local  current  densities  (as  seen  in  Fig.  5)  and 
creates  extra  stress  in  the  material  due  to  non-homogeneous 
volume  changes  in  the  electrodes,  leading  to  shorter  cell  lifetime. 

The  cell  voltage  as  a  function  of  the  amount  of  Li  in  the  positive 
electrode  is  presented  in  Fig.  7,  where  simulation  results  for  a  range 
of  discharge  currents  are  presented.  The  dashed  lines  represent  the 
reference  geometry,  while  the  solid  lines  represent  the  optimized 
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Fig.  6.  Concentration  distribution  in  the  electrodes  of  optimized  (a-c)  and  reference  (d-f)  3D-MB  geometries. 
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Fig.  7.  Discharge  curves  of  the  cells  with  optimized  and  non-optimized  3D-MB 
geometries  for  different  current  values.  Dashed  lines  represent  the  reference  geometry 
and  solid  lines  the  optimized  geometry. 


geometry.  It  is  seen  that  the  optimized  and  reference  geometry 
provide  very  similar  result  at  the  beginning  of  the  discharge. 
However,  when  the  lithiation  of  LixCo02  progress  beyond  x  =  0.75, 
the  performance  of  the  reference  cell  drops.  This  has  vital  influ¬ 
ences  on  the  battery  performance.  For  example,  the  discharge 
curves  for  the  optimized  geometry  at  4  A  m  2  and  the  non- 
optimized  at  2  A  m-2  discharge  current  density,  more  or  less 
overlap  for  x  >  0.85.  This  means  that  the  optimized  geometry  can 
provide  twice  as  high  current  for  the  same  voltage,  thus  signifi¬ 
cantly  increasing  the  power  of  the  battery. 

Generally,  the  discharge  curves  for  the  optimized  geometry 
remain  at  a  higher  voltage  for  a  longer  time  due  to  the  unchanged 
current  density  distribution  in  cell  (Fig.  5),  leading  to  uniform 
charging/discharging  of  the  electrode  materials  (Figs.  5  and  6).  In 
the  reference  geometry,  on  the  other  hand,  the  active  material  in 
the  plates  is  depleted  first,  followed  by  the  use  of  active  material  in 
the  electrode  bases.  As  a  result,  the  discharge  curves  of  the  opti¬ 
mized  and  reference  cell  coincide  during  the  first  half  of  the 
discharge  (x  <  0.75),  since  the  ion  transport  distances  are  similar. 
However,  the  ionic  transport  pathways  and  cell  resistance  increase 
when  the  charge  and  mass  transport  in  the  3D-cell  start  to  proceed 
mainly  through  electrode  bases  in  the  second  half  of  the  discharge 


cycle  (x  >  0.75),  resulting  in  a  cell  voltage  drop.  Similar  trends  can 
be  seen  also  for  other  discharge  curves. 

In  Fig.  8,  the  overpotential  and  resistance  of  the  cells  are  pre¬ 
sented.  The  dashed  line  represents  the  reference  (non-optimized) 
geometry  and  the  solid  line  represents  the  optimized  geometry. 
The  overpotential  (Fig.  8a)  follows  similar  characteristics  during  all 
discharge  currents,  obtaining  gradually  higher  values  when  the 
discharging  current  is  increased.  The  overpotential  in  the  optimized 
and  reference  geometries  almost  coincide  for  all  currents  for 
x  <  0.75  in  LixCo02.  However,  when  x  increases  above  0.75,  the 
effects  of  the  optimization  become  apparent.  The  overpotential 
required  to  operate  the  battery  increases  significantly  for  both  the 
reference  and  optimized  cell,  but  its  increase  is  significantly  lower 
for  the  optimized  geometry.  The  overpotential  tends  to  be  slightly 
lower  in  the  reference  geometry  at  the  beginning  of  discharge, 
when  x  <  0.65.  This  is  due  to  that  the  optimal  geometry  is  gener¬ 
ated  at  steady-state  conditions  which  are  not  yet  fulfilled  at  the 
beginning  of  discharge. 

The  cell  resistance  is  presented  in  Fig.  8b.  During  the  first  half  of 
the  discharge,  while  x  <  0.75  in  LixCoC>2,  the  resistance  is  almost 
independent  ofx  in  LixCoC>2  and  discharging  current  strength  for  both 
geometries.  Moreover,  the  resistances  are  almost  equal  for  both  cases. 
However,  when  the  discharge  reaches  the  second  half  of  the  cycle, 
beyond  x  =  0.75,  the  resistance  starts  to  increase  fast  for  both 
geometries,  corresponding  to  the  increase  in  overpotential.  The 
resistance  in  the  reference  geometry  is  almost  independent  of  the 
current  strength  in  the  second  half  of  the  discharge  cycle.  In  the 
optimized  geometry,  on  the  other  hand,  the  resistance  depends  more 
strongly  on  the  current,  and  is  higher  for  higher  current  values.  This 
can  be  explained  by  the  active  material  layer  thickness  at  the  tips  of 
the  electrodes,  which  leads  to  a  higher  concentration  overpotential  in 
the  electrolyte  in  the  electrode  pores  at  the  end  of  the  discharge  when 
higher  currents  are  used.  However,  this  effect  is  small  compared  to 
the  overall  decrease  of  the  overpotential  due  to  the  optimization. 

The  performance  increase  (defined  in  2.6)  due  to  geometry 
optimization  is  presented  in  Fig.  9.  Dashed  lines  represent  indi¬ 
vidual  performance  increases  for  different  discharging  currents  and 
solid  line  represents  average  performance  increase.  The  data  show 
the  amount  of  additional  energy  needed  in  the  reference  cell  (non- 
optimized  geometry)  as  compared  to  the  optimized  cell  to  maintain 
a  specific  current.  The  performance  of  the  optimized  cell  tends  to  be 
slightly  lower  than  the  performance  of  the  reference  cell  at  the 
beginning  of  the  discharge,  being  most  apparent  for  J  =  2  A  m  2. 
However,  for  other  current  values,  the  performance  increase  is 


Fig.  8.  Overpotential  (a)  and  normalized  cell 


i  (b)  for  different  discharge  currents  for  optimized  (solid  lines)  and  non-optimized  (dashed  lines)  geometries. 
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Fig.  9.  Performance  increase  in  the  optimized  cell  as  compared  to  the  reference  cell 
during  discharge. 


almost  the  same,  being  around  0.9.  The  average  performance 
increase  to  above  1  at  x  ~  0.63  and  remains  almost  constant  until 
x  =  0.75.  Within  this  range,  the  average  performances  of  the  opti¬ 
mized  and  non-optimized  cells  are  approximately  equal.  However, 
when  x  >  0.75,  the  performance  increases  fast.  The  average 
performance  increases  from  1  to  1.65,  indicating  that  the  energy 
needed  for  the  mass  transport  in  the  optimized  cell  is  65%  smaller 
than  in  the  reference  cell.  Moreover,  at  lower  current  densities 
(J=  2  A  m  2  and J  =  4  A  m  2),  the  maximal  performance  increase  is 
more  than  twice  as  high,  reaching  a  maximum  at  2.25  for 
J  =  4  A  m-2.  However,  when  the  current  density  increases,  the 
maximum  performance  increase  is  limited  to  1.4  for  J  =  10  A  m-2. 
This  performance  loss  at  higher  discharge  currents  for  the  opti¬ 
mized  geometry  is  caused  by  the  higher  concentration  over¬ 
potential  in  the  porous  electrodes,  seen  in  Fig.  8a. 

The  optimized  geometry  can  provide  as  much  as  twice  as  high 
current  as  the  non-optimized  one.  The  energy  dissipated  during  the 
discharge  in  the  cell  is  released  as  heat.  Thus,  up  to  two  times  less 
heat  is  released  due  to  the  optimization.  This  effect  increases  the 
safety  of  the  cells  by  decreasing  the  probability  of  overheating  and 
thermal  runaway.  The  lifetime  of  the  optimized  cells  will  also 
increase,  since  less  side  reactions  can  occur.  In  the  non-optimal 
geometry,  a  situation  can  occur  where  the  local  current  density  in 
certain  parts  of  the  electrode  (such  as  the  tips)  becomes  equal  or 
higher  than  the  limiting  current.  This  will  trigger  side  reactions 
which  can  destroy  both  electrode  and  electrolyte  material,  further 
reducing  cell  performance  and  ultimately  destroying  the  cell. 

It  is  important  to  address  the  question  whether  it  is  possible  or 
not  to  construct  this  kind  of  complex  non-uniformly  coated  elec¬ 
trode  structure  for  3D-MBs.  In  fact,  some  recent  studies  suggest 
that  non-uniform  coating  can  well  be  the  result  of  many  deposition 
techniques.  For  example,  Lethien  et  al.  [31  ]  recently  synthesized  an 
interdigitated  3D-MB  positive  electrode  of  silicone  nanopillars 
coated  with  layers  of  UPON  and  LiFePC>4,  resulting  in  a  geometry 
similar  to  the  one  presented  here:  the  pillars  tips  were  considerably 
thicker  than  the  rest  of  the  electrode.  Similar  results  were  also  seen 
in  studies  of  Lafont  et  al.  [32],  This  indicates  that  the  simulation 
results  presented  here  might  well  be  non-difficult  to  realize. 

4.  Conclusions 

In  this  paper,  a  variation  of  the  level  set  method  is  used  for 
optimization  of  a  novel  type  of  a  multiphysics  problem,  consisting 
of  coupled  diffusion  and  conductive  media  problems.  The  level  set 


method  is  used  to  optimize  the  electrode  shapes  in  a  3D-micro- 
battery  architecture.  The  optimized  geometry  was  thereafter  eval¬ 
uated  by  carrying  out  discharge  cycles  for  different  current  values. 
A  comparison  between  the  discharge  simulations  of  optimized  and 
reference  geometry  demonstrated  that  the  optimization  affects 
primarily  the  second  half  of  the  discharge,  when  x  >  0.75  in  Lix_ 
C0O2.  The  most  significant  effects  on  the  battery  voltage  develop¬ 
ment  occur  during  the  last  third  of  the  discharge  cycle,  from 
x  >  0.85  to  x  =  1.  The  optimized  geometry  demonstrates  consid¬ 
erably  better  results;  it  is  possible  to  reduce  the  energy  dissipated 
due  to  charge  and  mass  transport  up  to  2.25  times,  while  the 
average  performance  increase  was  up  to  1.65  times.  This  is  due  to 
the  fact  that  the  concentration  distribution  in  the  optimized  elec¬ 
trodes  demonstrates  a  uniform  lithiation/delithiation  change  over 
the  whole  electrode  during  discharge,  while  the  reference  cell 
electrodes  were  discharged  non-uniformly.  This  clearly  demon¬ 
strates  that  by  tailoring  the  material  deposition  onto  the  substrates 
when  constructing  3D-MBs,  higher  power  and  energy  density  can 
be  achieved. 
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Nomenclature 

a:  specific  area  of  an  electrode  (m-1 ) 
c:  concentration  (mol  m-3) 

D:  diffusion  coefficient  (m2  s_1) 

F:  Faraday  constant  (96,487  C  mor1) 

H,h,d.dc:  geometrical  parameters  of  the  cell  (Table  1 ) 
io:  exchange  current  density  (A  m-2) 

J:  ionic  flux  ^ 

r:  normalized  resistance  of  the  cell  (£2  m2) 

R:  universal  gas  constant  (8.314  J  (mol  K)_1) 

T:  absolute  temperature  (K) 

V:  volume  of  the  electrolyte  (m3) 

Voc:  open  circuit  voltage  (V) 
charge  number 
7*:  current  density  (A  m-2) 

t":  transference  number  of  Li+  ions  in  the  electrolyte 
If:  unit  normal  vector 


195  5(x):  Dirac  delta  function 

e:  porosity  of  the  electrodes 

£|S:  electrode-electrolyte  interface  thickness  during  optimization  (m) 
y:  surface  overpotential  in  the  battery  (V) 

X:  Lagrange  multiplier  (A  m“2) 
a:  conductivity  (S  m_1) 

<p:  electrical  potential  (V) 

£1:  geometrical  region  of  the  problem  (m) 
t:  pseudo  time,  representing  the  optimization  step 

Subscripts 
Li:  Li+  ions 
PF6:  PFg  ions 

0:  corresponding  to  constant  initial  value  or  boundary  condition 

i:  defined  when  used 
j:  defined  when  used 
l:  liquid  phase  (electrolyte) 
s:  solid  phase  (electrode  material) 
opt:  optimized  cell 

nonopt:  nonoptimized  (reference)  cell 


c p :  the  level-set  function 
a:  transfer  coefficient 
y:  reinitalization  parameter  (m  s_1) 


Superscripts 
n:  positive  electrode 
p:  negative  electrode 
*:  electrodes 


